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Integrability Tests for Nonlinear Evolution Equations 



1.1 Introduction 

During the last three decades, the study of integrabihty of nonhnear ordinary 
and partial differential equations (ODEs and PDEs) has been the topic of major 
research projects (see, e.g., [|], ^). This chapter presents a few symbolic algorithms 
to illustrate how computer algebra systems (CASs) can be effectively used in 
integrability investigations. We work with Mathematica [|42| , but our algorithms can 
be implemented in other languages. 

Among the many alternatives |Q for investigating the integrability of systems 
of PDEs with symbolic software, the search for conserved densities, generalized 
symmetries, and recursion operators is particularly appealing p^ . Indeed, it turns out 
that these quantities can be computed without the use of sophisticated mathematical 
tools. As a matter of fact, not much beyond differentiation and solving of linear systems 
is needed. As a result, our algorithms are easy to implement. In fairness, our algorithms 
are restricted to the computation of polynomial quantities of polynomial equations. 
Yet, this covers the majority of the cases treated in the literature. 

The algorithms in this chapter are based on a common principle: scaling (or 
dilation) invariance. Indeed, we observed that many known integrable systems are 
invariant under dilation symmetry, which is a special Lie point symmetry. The dilation 
symmetry can be computed by solving a linear system. Using dilation invariance, 
the plan is to first produce candidates for the polynomial densities, symmetries, and 
recursion operators in an efficient way. Once the candidate expressions are computed, 
their unknown constant coefhcients follow from solving a linear system. 

We focus our attention on explaining the strategy, at the cost of mathematical 
rigor and details, which can be found in jlj, |l^, ^ |2^. Rather than discussing 
the algorithms in general, we apply them to a few prototypical nonlinear evolution 
equations from nonlinear wave theory. Whenever appropriate, we address issues related 
to the implementation of the algorithms. For instance, we give explicit code for the 
Frcchet derivative, which is one of the key tools in our methods. 

Our package InvariantsSymmetries.m works for nonlinear evolution equations. 
Applied to a system with parameters, our package can determine the conditions 
on the parameters so that the system admits a sequence of conserved densities or 
generalized symmetries. Although we do not address it here, our package can also 
compute densities and symmetries of differential-difference equations (semi-discrete 
lattices). See O H for more information about that subject. 
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Due to memory constraints, our software can only compute a limited number of 
conserved densities and symmetries (half a dozen for systems; at best a dozen for scalar 
equations). To prove integrability, one must show that infinitely many independent 
densities or symmetries exist. Alternatively, one could construct the operator that 
connects the symmetries, and prove that it is a true recursion operator. Such proofs 
involve mathematical methods |@, ^ |3^, ^ that are beyond the scope of this article. 
Although it is not yet implemented, we also present an algorithm for the computation 
of recursion operators, based on the knowledge of a few conserved densities and 
symmetries. 

The computation of Lie point symmetries and generalized symmetries via 
prolongation techniques is purposely omitted. That topic and related software were 
covered extensively in [|l], . Space limitations also prevent the inclusion of the well- 
known Painleve test, which is a widely applied and successful integrability detector 
for nonlinear DDEs and PDEs. We refer to |Q for survey papers, books, and software 
related to Painleve analysis. 

This chapter is organized as follows. In Section 2, we discuss scaling symmetries 
of PDEs and show how to compute them. Section 3 deals with conservation laws. 
We give the definition and the steps of our algorithm, and show how to implement 
and apply the algorithm. We do the same for generalized symmetries in Section 4. 
An algorithm to determine the recursion operator is given in Section 5. The leading 
examples in Sections 2 through 5 are the Korteweg-de Vries (KdV) and Sawada-Kotera 
(SK) equations, and a system of nonlinear Schrodinger-type equations. For the latter, 
we derive the recursion operator in Section 6. In Sections 7 and 8, we discuss our 
Mathematica package InvariantsSymmetries.m and review similar software. We draw 
some conclusions in Section 9. 



1.2 Key Concept: Dilation Invariance 

Our algorithms are based on the following observation: if a system of nonlinear 
evolution equations is invariant under a dilation (scaling) symmetry, then its 
conservation laws, generalized symmetries, and the recursion operator have the same 
scaling properties as the system. This is at least true for the polynomial case. 

As leading example, we use the ubiquitous Korteweg-de Vries (KdV) equation ||l[, 

ut ^ Guuj; + U3^, (1.1) 

which describes water and plasma waves and lattice dynamics. Throughout this 
chapter, we will use the notations 

du d'^u d^u . , 



Equation (1.1) is invariant under the dilation (scaling) symmetry 

(t, X, u) (t/A^, x/A, \^u), (1.3) 



where A is an arbitrary parameter. Indeed, replacement of {t,x,u) according to 

J. 

dt 



(1.3) allows one to cancel a factor A^ in (1.1). Note that, e.g., ^ is replaced by 
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A'^^. Obviously, u corresponds to two derivatives in x, i.e., u ^ d^/dx^. Similarly, 
d/dtr^d^/dx^. 

We express all scalings in terms of Introducing weights, denoted by w, we could 
say that w{u) — 2 and w{T>t) = 3, if we set w(Dx) = 1- We used Dt and instead 
of w(^) and w{^) to cover cases where densities and symmetries depend explicitly 
on t and x (see |23| for examples) . 

The rank i? of a monomial is equal to the sum of all of its weights. Observe that 
(1.1) is uniform in rank since all the terms have rank R — 5, confirming that was 
a common factor. 

Computation of scaling symmetries. To compute the scaling symmetry of an 
equation, we compute the weights of all its terms, and require that the equation be 
uniform in rank. For (1.1), with w(T)x) = 1, this yields 

w(u) +u;(Dt) = 2u;(u) + 1 = w(u) +3. (1.4) 

The solution of this linear system is w(u) — 2, and w^Dt) — 3. 

As a second example, we consider a fifth-order PDE from soliton theory, 

Ut = ^u'^Ux + 5u3:U2x + 5uU3x + U^x, (1-5) 

due to Sawada and Kotera p8| . Scaling invariance requires that 

w{u) + w(Dt) = 3w{u) + 1 = 2w{u) + 3 = w{u) + 5. (1.6) 
Hence w{u) = 2 and w{Dt) = 5. 

Systems. Single PDEs like ( |l.l|) are a special case of 

Ut = F(U, U^,U22;, . . . ,U„a;), (1.7) 

where u and F are vector dynamical variables with n components. The number of 
components, the order m of the system, and its degree of nonlinearity are arbitrary. 

To determine the scaling symmetry, we require that each equation in (1.7) be 
uniform in rank, and solve the resulting linear system for the weights of all the 
variables. 

As an example, consider a vector nonlinear Schrodinger equation, 

Bt + (|BpB), + (Bo • B,)Bo + e x B,, = 0, (1.8) 

which occurs in plasma physic s ^. With Bq ~ {a,b) and B = (u,w) in the {y,z)- 
plane, and e along the x-axis, ( |l.8| ) can be written as 

Ut + [u(u^ + t;^) + /3u + 7f — Vx\ ^. — 0, 

vt+[v{u^ + v^) + eu + 5v + Ux]^^Q., (1.9) 

where [3 = a^,^ = 9 = ab, and S — are nonzero parameters. With reference to 



we call (1.8) or (1.9) the DMV equation. To start generally, we will consider the 



system ( |l.9| ) for arbitrary nonzero parameters (3,j,0 and S. 
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System (1.9) is not uniform in rank, unless we allow that the parameters [3 through 
5 have weights. Doing so, with ^(D^;) = 1, we obtain 

w{u) + w{V)t) = 3w(u) + 1 = w{il) + 2w{v) + 1 = w{u) + w(/3) + 1 

= w{v) + w(7) + 1 = w{v) + 2, (1.10) 

w{v)+w{T)t) = 2w{u) ^ w{v) + 1 = 'iw{v) + \ = w{u) + w{9) + \ 

= w(v) + w{5) + 1 = w{u) + 2. (1.11) 

Hence, 

u)(u) = w(w) = i, w{P)=w{-i)=w{e)=w{5) = l, w(Dt)=2. (1.12) 



Remark. For scaling-invariant equations like (1.1) and ( 1 . 5 ) , it suffices to consider the 
dilation symmetry on the space of independent and dependent variables. For systems 
like (1.9) that are inhomogeneous for scaling, we give weights to the parameters 
to circumvent the problem. For systems that lack scaling invariance and have no 
parameters, introducing one (or more) auxiliary parameter(s) with appropriate scaling 
provides a solution. 

The trick is to extend the action of the dilation symmetry to the space of 
independent and dependent variables, including the parameters. Doing so, our 
algorithms apply to a larger class of polynomial PDEs. The extra parameters are only 
used in the first step of the algorithms: that is, in producing the candidate densities 
and generalized symmetries. Beyond that first step, parameters are no longer treated 
as dependent variables! Details and examples are given in ^ p3| . 



1.3 Conservation Laws 



Definition. A conservation law for (1.7), 

Dtp + D,J-0, (1.13) 
connects the conserved density p and the associated flux J. As usual, D^ and D^, are 



total derivatives, and (1.13) holds for all solutions of (1.7). Hence, density-flux pairs 



only depend on u, u^, etc., not on Uf . With a few exceptions, densities and fluxes do 
not explicitly depend on t and x. 

For the scalar case, Ut = F, the computations are carried out as follows: 



dt ^-^ dUkr 
k=0 



(1.14) 



where n is the order of p. Upon replacement of ut,Uxt, etc. from ut — F, one gets 



DtP=^+p'(^^)[F], 



(1.15) 



where p'{u)[F] is the Frechet derivative of p in the direction of F (see Section 1.4). 
Furthermore, 



D,J: 



dJ 



E 



dJ 



dx duk. 



U(fc+l)a;i 



(1.16) 
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where m is the order of J. Integrating both terms in ( |1.13D with respect to x yields 



/ + 00 
pdx = -J\tZ = Q, (1.17) 
-oo 



provided that J vanishes at infinity. In that case, 

/ + 00 
p da; = constant in time. (1-18) 
-oo 

So, P is the true conserved quantity. For ODEs, the quantities P are cahed constants 
of motion. 

Examples. The first three (of infinitely many) independent conservation laws [0, |29| 
for (|]|) are 

Dt(u)- 0,(3^2+^,2^)^0, (1.19) 
Dt(u2) - D,(4u3 -ul + 2uu2^) - 0, (1.20) 

Dt {^^~ ^^^^ ~ (^"^~ 6mM^+ 'iu'^U2x+ ^^2,- U,U3,^ = 0. (1-21) 

The first two conservation laws correspond to conservation of momentum and energy. 
Note that the above conservation laws are indeed invariant under ( |l.3D . The terms in 
the conservation laws have ranks 5, 7, and 9. The densities 

p(i)=K, p^^^=u\ and p(3)=i.3_1^2 (^ 22) 

have ranks 2, 4 and 6, respectively. The associated fluxes have ranks 4, 6 and 8. 
Equation (1.1) also has a density- flux pair that depends explicitly on t and x : 

p ^ tu^ + ^xu, (1.23) 

J = t{4:u^ + 2uu2x-ul)+x(^u^ + ^U2x^ ~^Ux. (1.24) 

p has rank 1, J has rank 3, since w(t) ~ —3 and w{x) = —1. To accommodate this 
case, w e us ed the total derivative notation Dt and D^,, in ( 1.13| ). 

For (|l.5|) , the first two (of infinitely many) conserved densities ||l^ are 

= u and ^ 1 3 _ ^2 (^ 25) 



We will use them in the construction of the recursion operator for (1.5) in Section 5. 

We now describe how to compute polynomial densities that are (explicitly) 
independent of t and x. We refer to |l^ for the general algorithm covering systems 
as well as {x, i)-dependent densities. 

Algorithm for Polynomial Conserved Densities 
Step 1: Determine the form of the density 
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Select the rank R of p, say, R — 6. Make a list C of all monomials in the components 
of u and their x-derivatives that have rank R. Remove from C all monomials where 
the power of the highest derivative is 1. This is done to remove terms in p that are 
in Image (Dx), and therefore belong to the flux J. After all, densities are equivalent if 
they only differ by terms that are total derivatives with respect to x. 

Make a linear combination with constant coefficients Ci of the monomials that 
eventually remain in the list C. 



For (LI), C = {u ,u^,uu2x,U4x}- Next, u^x and uu2x are removed. Obviously, 
U4x = ^xUsx, and uu2x — ^D^w^ — u^. So, uu2x and only differ by a total x- 
derivative. From C = {u^, u^}, one constructs p = ciu^ + C2it^, which has rank R = 6. 

Step 2: Determine the unknown coefficients 

Substitute p into the conservation law ( |1.13| ), and compute Dtp via ( 1.14 ). Use the 



PDE system to eliminate all t-derivatives of u, and require the resulting expression E 
to be a total x-derivative. 

To avoid integration by parts, apply the Euler operator (also called the variational 
derivative) Q 



dUkx 

d ^ f d \ ^.f 



to E of order m. liLu{E) = immediately, then i? is a total x-derivative. If L„(i?) ^ 0, 
then the remaining expression must vanish identically. This yields a linear system for 
the constants c^. Solve the system. Carrying out these operations for (^), one gets 

Ci = 1,C2 = -\. 



Remark. With (1.15), the system for the Ci follows from Lu{p' {u)[F]) = by 



equating to zero the coefficients of monomials in u and their x-dcrivatives. 
Implementation in Mathematica 

In Mathematica^ D is the total derivative operator and the variational derivative (Euler 
operator) can be found in the Standard Add-on Package Calculus'VariationalMethods' . 



For instance, returning to p^ ' in ( 1.22 ), with (LI) and ( 1.14 ), one computes 



E = Dtp(3) ^ p(3)'(u)[Mt] = (3u' - UxB)ut 

— ISu'^Ux — 6u'^ ~ 6uUxU2x + 'iu'^U^x — 'U'xU4x. (1-27) 

Application of the variational derivative, VariationalD [E,u[x,t] ,-[x,t}] gives zero. 
That means that i? is a total x-derivative of a polynomial J^'^^ Integration of 
E = -BxJ^^^ gives 

j(3) ^ „(^y4 _ gy^2 ^ 3^2^^^ ^ 1^2^ _ y^yg^). (1.28) 

Example. With our package InvariantsSymmetries.m, we searched for conserved 



densities of ( |1.9D . Obviously, (1.9) is a conservation law; thus, p*-^-* = u and 
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without conditions on the parameters. Additional conserved densities only exist if 
^ — 0. The first few are: 



(3) = „2 



p ' = 



V 

(4) ^ L{u^ + y^f + {j3- 5)u'^ + 2euv + 2vu^, 



+ - 5){u* - V*) + 'iu^VUx + V^U^c, 



(1.29) 
(1.30) 

(1.31) 



and 



32 



in' 



- V"")^ + + f')(?4 + Vl) + \{UU^ + V^x? 



1 



'2„,2 



- '5)^^« - le^uv + 'Uu^v^ - le\u, 

O L L L 

3 1 11 

Via integration by parts, p^^') is equivalent to 



(1.32) 



(1.33) 



Indeed, p^^) = pi^) + D^i 



Judged from ( 1.2£ )-( 1.32 ), the complexity of the expressions dramatically increases 
as the rank increases. The fact that we were able to compute 6 independent densities 
for (1.9) is an indicator that the system presumably is completely integrable, as was 
later proved in . 

Only the conserved densities p*-^-* through p*-"^) Tf^m used in the construction of 
the recursion operator for (l.£) in Section 6. 



1.4 Generalized Symmetries 



Definition. A vector function G(a;, t, u, u^,, U2x, • • •) is called a symmetry of (1.7) if 
and only if it leaves (^^) invariant under the replacement u — » u + eG within order 
e. Hence, 

Df(u + eG) =F(u + eG) (1.34) 



must hold up to order e on any solution of (1.7). Consequently, G must satisfy the 
linearized equation ^ 

DtG = F'(u)[G], (1.35) 
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where F' is the Frechet derivative of F, i.e 

d_ 
de 



F'(u)[G] --F(u + eG)|,=o- (1-36) 



In (1.34) and (1.36), we infer that u is replaced by u + eG, and Unx by Unx + eD"G. 



As usual, Dt and D^, are total derivatives, and G = (Gi, G2, . . . , G„) if the system 



(1.7) has n components. 

Once higher-order symmetries have been found, these vector fields can be used to 
obtain fundamental information about the integrability of the equation. In many cases, 
conserved quantities, Hamiltonian structures, and recursion operators follow readily 
from the knowledge of generalized symmetries W . 



Examples. The first three (of infinitely many) symmetries |3^] of (1.1) are 

G(1) = Ux, G(2) =6uUx+U3x, 

G^^^ = 30u^Ux + 20uxU2x + Wuu3x + U5x- (1-37) 

All the terms in these symmetries have rank 3, 5 and 7, respectively. 

With higher-order symmetries one can generate new integrable PDEs. For example, 
Ut = G^'^^ is the Lax equation in the completely integrable KdV hierarchy jl]. 



Note that (1.1) also admits symmetries that explicitly depend on t and x. Indeed, 
the symmetries 

= 1 + 6tux and G^^^ ^4u + 2xux + 36tuux + 6tu3x (1.38) 

are of rank and 2. They linearly (and explicitly) depend on t and x. 

The algorithm presented in this paper can easily be extended to cover this type of 
symmetries (see |1J, \l% for details) . 



The situation for (1.5), which also has infinitely many polynomial symmetries, is 



more complicated. The symmetries of (1.5) originate from two distinct "seeds": 



G'^^^ — Ux and G*-^-* = 5u^Ux + 5uxU2x + ^uu^x + u^x- (1.39) 

We have also computed symmetries of higher rank, but we do not show them here. A 
detailed computer-aided study showed that the symmetries G^^'~^^ with rank 6i — 3 
come from the seed G'^^-', whereas G*^^*^ with rank 6i -I- 1 originate from G'^^-', where 
i — 1,2, . . . (sec Q for details). 

For systems of type (^]^), the symmetry G is a vector with n components. Our 



computer search with InvariantsSymmetries.m revealed that (1.9) is invariant under 
the transformation (u,v) — ^ {v,—u), which is a Lie point symmetry, provided the 
conditions (3 = 5 and 7 = —0 hold. However, these conditions do not lead to a 
hierarchy of integrable equations. We therefore continued our search with arbitrary 
nonzero parameters (3 through 5. 

The first two symmetries of (0) are G^^) = (gJ^',G^^') and G^^) = (Gf\G^^^), 
where 

G\ = Ux, 

G^'^ = vx, (1.40) 

G^-' = iP - 5)Ux + 'iu^Ux+v'^Ux'^^Vx + 2uVVx-V2x, 

G^2^ = 9ux + 2uvUx + u^Vx + iv'^Vx + U2x- (1-41) 
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Note that the sum of symmetries is still a symmetry. Remembering (ut, vt) — (-Fi, F2) 

'2 



IJl), we then have G^^ + 5G^^^ = -Fi and G^^^ 



The next symmetry, G^'^) = {G\ , G2 ), only exists if 7 = It is 

15 3 
Gp) = 3(/3 - 5Wux + -u\x + 60UVUX + 9uWux + -v\x + ie{u' + v^)vx 

-(-6(m'^w + uv^)vx — 3(u^+ v'^)v2x — QuUxVx — 6vv'^ — U3x, (1-42) 
G^^ = 39{u' + v^)ux + 6{u' + v^)uvux + 6uul + ^u^v^ + WuvVx 

15 

■A^U^V^V^ - 3(/3 - 6)v'^Vx + -l^'^'^'^x + QvUxVx + 3(u^+ W^)u2a; - W3a;- (1-43) 

The algorithm for symmetries is similar to the one for conserved densities. The only 
difference is that monomials that differ by a total x-derivative are no longer removed 
from the list C. 

Algorithm for Polynomial Generalized Symmetries 
Step 1: Determine the form of the symmetry 

Select the rank R of the symmetry. Make a list C of all monomials involving u 
and its a;-derivatives of rank R. To obtain the form of the symmetry, make a linear 



combination of these monomials with constant coefficients Ci. For example, for (1.1), 
G = ci v?Ux + C2 UxU2x + C3 uu^x + C4 u^x IS the form of the generalized symmetry of 
rank R = l. 

Step 2: Determine the unknown coefficients 

Compute DfG. Use the PDE system to remove all i-derivatives. Equate the result 
to the Frechet derivative F'(u)[G]. Treat the different monomial terms in u and its 
^^derivatives as independent to get the linear system for the c^. Solve that system. For 



(1.1), one obtains the symmetry of rank 7 : 

G = ZQu^Ux + 20uxU2x + Wuu^x + u^x- (1-44) 



Symmetries of lower or higher rank are computed similarly. See [g_4| |lj, ^ for details 
about the algorithm and its implementation. 

Remark. Starting with a conserved density p, the symmetries for a Hamiltonian 



system can be obtained from 'Dx{Lu{p)), where is defined in (1.26). See, for 
example, for a study of the connection between densities and symmetries for 
Lagrangian and non-Lagrangian systems. 

Implementation in Mathematica 

The key tool to compute symmetries is the Frechet derivative, which is implemented 
as follows: 

frechet [f uncF_List , f uncU_List , indVars_List , f uncG_List] : = 
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Module [{eps , resultlist={} , i} , 
Do [resultlist=Appeiid [resultlist , 
Expand [(D [Part [funcF.i] /. 

{(f_/; MemberQ [funcU.f] ) [Sequence @@ indVars]:> 

f [Sequence 00 indVars] +eps*funcG [ [Flatten [Position [funcU, f] ] [[1]]]] , 
Derivative [k_ , 0] [f_] [Sequence @@ indVars] :> 

D[f [Sequence @@ indVars]+eps*funcG[[Flatten[Position[funcU,f]] [[1]]]], 
{indVars [ [1] ] , k}] } , eps] / . eps : > 0) ] ] , 
{i , Length [f uncF] }] ; 
Return [resultlist] ] ; 

To compute the Frechet derivative of — ^u^ in the direction of k{x,t), type 
frechet [{u[x,t] -3-(l/2)*D[u[x,t] ,x] '2},{u>,{x,t},{k[x,t] }] ; 
This gives the answer 3u^k — Uxkx- 



1.5 Recursion Operators for Scalar Equations 

Definition. A recursion operator is a linear operator $ on the space of differential 
functions with the property that whenever G is a symmetry of ( p_.7| ), so is G with 
G = *G. 

The equation for the recursion operator |3^, |3^ is 

Dt* + F'{u)] = -g^ + *'[F] + * o F'{u) - F'{u) o * = 0, (1.45) 

where [ , ] means commutator, o indicates for composition, and the variational 
derivative of the operator $ is defined in, e.g., 0|. 

For 71-component systems like (^3), the symmetries G are vectors with n 
components, and the recursion operator $ is an n x n matrix. 



Examples. The recursion operator for (LI) is 



$Kdv = + 4u + 2ux'D-^ = 



2u + 2DuD" 



(1.46) 



where, for simplicity of notation, D — and D^^ = D~^ . Here, $Kdv G*^^^ = G'-^-' 
and <&Kdv G^^^ = G^"^^ for the symmetries listed in ( 1.37 ). The recursion operator 
<i>Kdv also connects the (x, t)-dependent symmetries ( 1.38 ), i.e., ^KdyG*-^^ — 
but ^KdvG'^' is no longer local. 

Since w(D~^) = — w(D) = —1, the three terms in ( 1.46| ) have rank R — 2. The 
recursion operator ( 1.46| ) is uniform in rank. Clearly, the rank of <I>Kdv is the difference 
in rank between consecutive symmetries in (1.37). 

In view of the symmetries ( I.39|) , the recursion operator for (1.5) must have rank 
6. Indeed, the recursion operator ||9|, |l^ has rank 6 : 



$SK =D^ 



2uD* + 2DuD^ 



3uDmD 



SuB^u - 2DuDm - 2u^ 



5DmD^mD"^ 



2u^B) 



(1.47) 
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which can also be written as 



$SK = + 3ul)^ - 3DuB^ 



llD^uD^ - lOD^uD + 5D^u 



+12u2d2 - 19uDuD + SuD^u + 8DuDu + 



(1.48) 



with G(2) in (jTsI). 

Our algorithm for the computation of polynomial recursion operators is based on 
the following observations. 



Key observations. All terms in ( 1.46 ) and ( 1.48 ) are monomials in D, D~^,u, and 
Ux- Depending on the form of the recursion operator, U2x,u^x, etc. can also appear, 
as is the case in ( 1.47 ). 

Recursion operators split naturally in $ = $0 + $1, where $0 is a differential 
operator (without terms), and $1 is an integral operator (with terms). 

Furthermore, application of <I> to any symmetry should not leave any integrals 
unresolved, since all symmetries are polynomial (see [Q). This is where the connection 
between conserved densities and symmetries comes into play. 

For instance, for ( |l.l| ) it is clear that T)~^ {QuUx + ws^^) — + U2x is polynomial. 
Similarly, for ( |l.5[ ) , using ( |l.4g ) , we have 



D ^{bu^Ux + 5UxU2x + 5uU3x + U^x) = -U^ + 5uU2x + U4x 



(1.49) 



and 



D ^{u^ - 2UxD)i5u'^Ux + 5UxU2x + buUSx + Ur,x) 



The first two conserved densities of (1.5) are p*^^^ = 



u and p^^^ — ^u"^ — u 



(1.50) 
l. Thus, with 



( 1.13 ), we get DfM — ut — — D^; j(^) and 

Dt - ul^ = p'{u)[ut] = {u^ - 2uxB)ut = -D,J(2). (1.51) 



So, the factor {u^ — 2uxTy) in ( 1.48 ) comes from p^'^\ and D^^[(u^ — 2uxB)ut\ will be 
polynomial, namely ~J^^\ 

A similar situation happens for (1.1), where p'-^-' = u, p'^- 
Then, with ( |l.l3D and (|lT5| ), 



i^, and p''^) 



2^x- 



Dtp(i) = D 



tu = Ut 



2uut 



'J:>xJ^'^\ and 



Dtp(3) = Di(u3 - \l) = p(3)'(u)K] = {iu" - UxDx)ut = -D, J(3), (1.52) 

for polynomial J*^*', i — 1,2,3. Thus, application of D^^, or D^^w, or D^^(3u^ — u^D) 
to QuUx + u^x leads to a polynomial result. However, as will be shown below, the 
terms involving D~^u and D^^(3w^ — w^D) are not needed in the construction of the 
recursion operator (1.46|). 
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Since our algorithm for recursion operators is the most elaborate, we give the steps 
that lead to (Oel) and (p 



1. 



We consider the scalar case first. Systems are dealt with 



in Section 6. 



Algorithm for Polynomial Recursion Operators 
Step 1: Construct the form of the recursion operator 
(i) Determine the rank of the operator 

Compute the rank R of the operator based on the known ranks of consecutive 
symmetries. For example, from (|l.37 ), we compute 

R = rank $ ^ rank G^^^ - rank G^^^ = rank G^^^ - rank G^^^ = 2. 



(1.53) 



Obviously, rank $o = rank $1 ~ rank ^ — R. 



(ii) Determine the pieces of the operator $0 

Make a list C of all permutations of D-'m'^, with j and k nonnegative integers, that 
have the r ank R. For (pl), 

/: = {D^m}. (1.54) 



(iii) Determine the pieces of the operator $1 

It can be shown [|[ ^ that 

*i=EEG^'^D-V'=)'(«), (1.55) 

J k 

where the symmetries G'--'^ are combined with and p^*^-* (u) in such a way that 
every term is exactly of rank <I>i — R. That is, the indices j and k are taken so that 
rank (G*^-'^) + rank {p^''^ (u)) — 1 = i? for every term in ( 1.55 ). 

Using the densities and symmetries, make a list M. of the pieces involving D^^. 

For (LI), from Table [| it should be clear that can only be sandwiched between 
Ux and 1. Any other combination would exceed rank 2. Hence, 

X = {u^D-i}. (1.56) 

(iv) Build the operator $ 

Next, produce TZ = C[JA4, which has the building blocks of the recursion operator. 
To get $, linearly combine the pieces in TZ with constant coefficients c^. For (p~l|), we 
obtain 

7^= {d2,u,m,D-1}. (1.57) 

Thus, 

*Kdv = ci + C2 u + C3 u^D-^ (1.58) 



We now repeat steps (i)-(iv) for the SK equation (1.5). 
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Table 1 Building blocks of $i for the KdV equation. 



Rank 


bymmctry G^-" 


Density ' 


' (U) 



1 






1 


2 




u 


u 


3 


Ux 






4 








5 








6 









(i) Using the symmetries ( 1.39 ), we get 

rank $ = rank G^^^ - rank G^^^ = rank G^^^ - rank G^^^ 6. (1.59) 

(ii) The operator $0 will be built from 

C = {D^uD^,DuD3,D2mD2,D3uD,D^w,u2d2,uDmD, 

MD2M,DM^D,DMD^i,D2^i^^t^}. (1.60) 



Table 2 Building blocks of $1 for the SK equation. 



Rank 


Symmetry G*--'-' 


Density p'-'^-' 


p^^y{u) 



1 
2 
3 
4 
5 
6 
7 


hu^Ux + ^UxU2x + ^UU^x + U5x 


u 

_ ul 


1 

V? — 2uxO 



(iii) From Table which list the building blocks for <i>i for (1.25), we obtain 

M = {uj;D"^(m^ - 2m^D), {bu^Ux + hUxU2x + ^uu^x + U5x)D"^}- (1-61) 
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All other combinations of the form G^-'-'D^^p'^'^'' (u) exceed rank 6. 
(iv) Combining the monomials from TZ ~ C[J A4, we get 

$SK = ciD^ + C2uD^ + csDuD^ + aD^uD^ + csD^uD 

+060"*^ + ctm^D^ + csmDuD + CguD^M + CioDu^D 

+ Ci5(5u^U2; + 5UxU2x + 5uU3x + M5a;)D""^- (1-62) 



Step 2: Determine the unknown coefficients 

To determine the coefficients Ci, require that 

^Q(k) ^ Q{k+s) ^ A: = 1,2, 3,. 



(1.63) 



where s is the number of seeds. In practice, it suffices to use k — 1 and 2 in ( 1.63 ) to 
fix all coefhcients c;. Solve the resulting linear system(s) for the unknown 

For dn]), $KdvG^^^ = G^^'> with $Kdv in ( [L5^ ) and its symmetries ( [OTI ), we 
obtain 

5 = {ci - 1 = 0, 18ci + C3 - 20 = 0, 6ci + C2 - 10 = 0, 2c2 + cg - 10 = 0}. 



The solution is ci = 1, C2 = 4, and C3 = 2. Substituting it into (1.5S), we get 

$Kdv = D2 + 4u + 2u^D-\ (1.64) 

The ex plic it computation on page 260 in (t) shows that ( 1.64 ) satisfies ( 1.45 ). 
For ( L5 ) , we express that 



ci>SKGW=G(3) and ^skG^^) = G^^), 



(1.65) 



with <i>sK in (1.62) and the symmetries (1.39). Then we solve for the constants c^. This 
eventually yields 



$SK = + 3mD* - 3DuD^ + llD^iiD^ - IOD^mD + 5D^u 
+12u2d2 - 19uDuD + S-itD^-it + 8DmDu + 4u^ 



(1.66) 



with G^^-* in (1.39). A lengthy computation shows that this recursion operator satisfies 

( psD - 

After integration by parts, (1.47) or, equivalently, ( 1.66| ) can also be written |Q as 

$SK = + 6mD'' + 9m:^,D''^ + gu^D^ + llu2xO'^ + lOua^^D + 21uUxD 

+4u^+ 16uM2x+ Gul+ 5u4x+ UxB-'^{u'^+ 2u2x)+ G^^^D^^ (1.67) 



1.6 Recursion Operators for Systems 

We show how to construct the recursion operators for systems (^]^) with n components. 
The symmetry G has n components and the recursion operator €> is a n x n matrix. 
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We used Mathematica interactively to compute the recursion operator for ( |l.9| ) with 
7 = 0. The recursion operator has the form 



(1.68) 



$11 $12 
$21 *22 

^From G*^^' = ^G^"'^', the rank of the entries $y is determined by 

rank G'l^'' = rank $ii + rank G^^^ = rank $12 + rank Gj^"* , 

rank G^^^ = rank $21 + rank g[^^ = rank $22 + rank G^^^ . (1.69) 

For ( |l.9| ), the difference in rank between consecutive symmetries is 1, so rank $ij = 
1, i,j = l,2. 

(i) As for th e sca lar case, we first construct the differential operator $o- In view of 
the weights ( 1.12 ), 

C,j = {^t^ uv, D, /?, S, 6}. (1.70) 

Hence, 

^ _/ CiU^+C2W^+C3WW + C4D + C5 Ceti^+CyW^+CgMW + CgD + Cio A (111) 

° \^CiiU^+Ci2W^+Ci3UU + Ci4D + Ci5 CiGu'^+Ci7v'^+CisUV + CigB + C2oJ ' ^ ' ' 

where C5, Cio, C15, and C20 will be linear in /3, 6, and 9. 

(ii) Using the conserved densities p^^^ — u, p'-^^ — v, and p^^^ — v?' + , we have 

DtpW = Dtu=|^^^i + |^F2 = (l,0)-(ui,«0 = -D.^''\ (1-72) 
Dip(2) = Dii; = |^^^i + |^^^2-(0,l)-(7/t,i;0 = -D.>/^'\ (1-73) 

Dtp(3) ^ Dt(u2 + «2) = ^_ Lf^ + ^- IF2 

ou ov 

= 2{u,v) ■ {uuvt) = -'D^J^''\ (1.74) 

where the dot (•) refers to the standard inner product of vectors. Therefore, introducing 
the symmetry {ux,Vx) on the left of gives 

M = [{ux,vx)^ QTi~\u,v)Y (1-75) 
where © stands for the tensor product of matrices and T for transpose. So, 



*i = C2i{ux,Vx) 0D ^(u,w) 

C2lUxT)^^U C2lUx'D^^V 
C2lVxD^^U C2lVxD^^V 



(1.76) 



Note that {ux,Vx)'^ D ^(1, 0) and {ux,Vx)^ D ^(0, 1) are of rank ^. They cannot 
be used in $1, where all pieces must have rank 1. 
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To determine the unknown constants in $ = $o 



k 

1=1 



- $1, we use 
k 



1,2,. 



(1.77) 



where aki are unknown coefRcients (which can be zero). In contrast to the examples 
in the previous sections, the aki play a role when dealing with systems with weighted 
parameters like (1.9). 

It suffices to take k — 1 and 2 in (1.77) to fix all coefficients Ci, and the extra 
unknowns an, 0:21 and a22- By solving a linear system, we obtain an = 0,a2i = 0^ 
and a22 = f3 — S, and the values for the coefficients Ci through C21. The recursion 
operator then follows readily: 



/3 - (5 + 2u2 - 
9 + 2uv + D 



2^;^D"l^t 



* = 

The recursion operator for the case j ^ 9 



2uv 
2^2 



D + 2u^D-iw 



(1.78) 

was computed analytically in WOl. 



1.7 About the Integrability Package InvariantsSymmetries.m 

We briefly describe our package InvariantsSymmetries.m, which automatically 
performs the computation of conservation laws (invariants) and symmetries based 
on the algorithms in Sections 3 and 4. 

Users must have access to Mathematica 3.0 |Q. The files for our package are 
available from MathSource ||l8|. The package includes instructions for installation, 
on-line help, documentation, and built-in examples. 

After proper installation, it is advisable to first run our notebook (called Examples), 
which is accessible through the browser as part of the Add-on Package Integrability. 
The interactive examples in the notebook will help familiarize the user with the syntax 
of our functions (see also [|l8j). 

To use the package as part of a new notebook, start Mathematica and type 
In[l]:= <<Iiitegrability ' to read in the package. You will get the following 
message: 

Loading init.m for Integrability from AddOns . 

The key functions for conservation laws and symmetries of PDEs are PDEInvariants 
and PDESymmetries. These functions take the following arguments: the equations in 
the system, the dependent and independent variables, and the range for the rank. 



For example (1.9), the first two lines below define the system (|l.9|). The third line 



(1.42) 



produces the densities (1.29)-(1.31). The fourth line gives the symmetries (1.40) 



In[2]:= pdel := D[u[x,t] ,t]+D[u[x,t]*(u[x,t] ~2+v[x,t] "2) + 
beta*u [x , t] +gcimma*v [x , t] -D [v [x , t] , x] , x] == ; 

In[3]:= pde2 := D[v[x,t] ,t]+D[v[x,t]*(u[x,t] "2+v[x,t] '2) + 
theta*u[x,t]+delta*v[x,t]+D[u[x,t] ,x] ,x] == 0; 



whugchap 9/2/2008 04:31— PAGE PROOFS for John Wiley & Sons Ltd (31x47jw.sty v5.0, 15th April 1997) 



INTEGRABILITY TESTS 



17 



In[4] := PDEInvariants[{pdel,pde2}, {u,v}, {x.t}, {1,3}]; 

In[5] := PDESymmetries[{pdel,pde2>, {u,v}, {x,t>, {3/2,5/2}]; 

Help about the functions and their options is available on-line. For instance, type 
In [6] := ??PDEInvariarLts to obtain the function description. Part of it reads: 

PDEInvariants [eqn, u, {x,t}, R, opts] finds the invariaint with raink R 

of a partial differential equation for the function u. 

PDEInvariants [{eqnl,eqn2, .. .}, -[ul,u2, . . .}, {x,t}, {Rmin , Rmax} , opts] 

finds the invarieints with rank Rmin through Rmax. 

X is understood as the space variable and t as the time variable. 

Typing In [7] := ??PDESymmetries produces descriptions like: 

PDESymmetries [eqn, u, {x,t}, R, opts] 

finds the symmetry with rank R of a partial differential equation 
for the function u. 

Information about the options of PDESymmetries is obtained by typing 
In [8] := ??WeightedParameters. It returns: 

WeightedParsuneters is an option that determines the parameters with 
weight. If WeightedParamieters -> {pl,p2,...}, then pi, p2, ... are 
considered as constant parameters with weight. 
The default is WeightedParameters -> {>. 

The option WeightedParsuneters is useful when working with systems that lack 
uniformity in rank. In such cases, our software tries to resolve the problem by itself 
and prints appropriate messages. When unsuccessful, the program will suggest the 
use of the WeightedParameters option. Therefore, the option WeightedParameters 
should not be used until it is explicitly recommended by the software. 

Rules for the weights of variables can be entered via the option WeightRules: 

WeightRules is an option that determines the rules for weights of 
the variables. If WeightRules -> {Weight [u] -> val,...}, then scaling 
properties are determined under these rules. There is a built in 
checking mechanism to see if the given rules cause inconsistency. 

For PDEs, the MaxExplicitDependency option allows one to compute conserved 
densities or symmetries that explicitly depend on the independent variables: 

MaxExplicitDependency is an option for finding the invariant and 
generalized symmetries of PDEs and DDEs. 

If MaxExplicitDependency -> Max_Integer, then the program allows for 
explicit dependency of independent variables of maximum degree Max. 
The default is MaxExplicitDependency -> 0. 

1.8 Software Review 

In this section, we briefly review software for the computation of conservation laws, 
higher-order symmetries and recursion operators. In Table ??, we give a summary and 
contact information. 
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Higher-order symmetries can be computed with prolongation methods, and 
numerous software packages are available that can aid in the tedious computations 
inherent to such methods. A 50 page survey of software for Lie symmetry 
computations, including generalized symmetries, can be found in pl| ], and a short 
update in ||2^] . We will not repeat these software reviews here. A survey of packages for 
conservation laws was first given in ||l5|] . However, to keep this chapter self-contained, 
we present a summary of that survey. 



Based on dilation invariance, Ito's programs in REDUCE (see [g5|, |26| |27|) compute 
polynomial higher-order symmetries and conserved densities for systems of evolution 
equations that are uniform in rank (no weighted parameters can be introduced). Ito's 
latest program, called SYMCD, cannot be used to compute symmetries and densities 
that depend explicitly on the independent variables t and x, nor can it handle systems 
with parameters. More details are given in [l^ . 

In I ^, ^ , Fuchssteiner et al. present algorithms to compute higher-order symmetries 
of evolution equations. Their algorithm in |^ is based on Lie- algebraic techniques 
and uses commutator algebra on the Lie algebra of vector fields. Their approach is 
different from the usual prolongation method in that no determining equations are 
solved. Instead, all necessary generators of the finitely generated Virasoro algebra 
are computed from one given element by direct Lie-algebraic methods. Their code 
is available in MuPAD. In pOl, Fuchssteiner et al. give code to verify that recursion 
operators are hereditary. In [pj, it is shown how to compute mastersymmetries from 
which the recursion operators can be retrieved. 

The REDUCE program FS for "formal symmetries" was written by Gerdt and 
Zharkov (see also 1 1 1, 12j ) . FS computes higher-order symmetries and conservation 
laws of polynomial type. The algorithm requires that the evolution equations be of 
order two or higher in the spatial variable. However, this approach does not require that 
the evolution equations be uniform in rank. With FS, one cannot compute symmetries 
that depend explicitly on the independent variables t and x. Applied to equations 
with parameters, FS computes the conditions on the parameters using the symmetry 
approach. 

The PC package DELiA, written in Turbo Pascal by Bocharov ||] and co-workers, 
is a commercial computer algebra system for investigating differential equations using 
Lie's approach. The program deals with higher-order symmetries, conservation laws, 
integrability and equivalence problems. It has a special routine for systems of evolution 
equations. The program requires the presence of second- or higher-order spatial 
derivative terms in all equations. For systems with parameters, DELiA does not 
automatically compute the densities and symmetries corresponding to the (necessary) 
conditions on the parameters. One has to use DELiA'' s integrability test first, to 
determine the conditions. Once the parameters are fixed, one can compute the densities 
and symmetries. 

Sanders and Wang have Maple and FORM code for the computation of symmetries 
in the scalar case, allowing zero and negative weights p6| , ^ and nonpolynoniial 
equations and symmetries. This code relies on the Maple package dijfalg Q| to do 
the reductions of solutions of ODEs (PDEs). See |^ for theoretical foundations of 
the computation of conservation laws, and for the use of their algorithms in the 
integrability classification of KdV-type higher order PDEs. 

Wolf et al. H have three packages, called CONLAW 1/2/3, in REDUCE for the 
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computation of conservation laws. There is no limitation on the number of independent 
variables. The approach uses Wolf's program CRACK for solving overdetermined 
systems of PDEs (see [^l], ^). Wolf's algorithm is particularly efficient for showing 
the non-existence of conservation laws of high order. In contrast to our program, it 
also allows one to compute nonpolynomial conservation laws. 

Hickman ||2^ at the University of Canterbury, Christchurch, New Zealand 
has implemented a slight variation of our algorithm for conserved densities in 
Maple. Instead of computing the differential monomials in the density by repeated 
differentiation, Hickman uses a tree structure combining the appropriately weighted 
building blocks. 

1.9 Conclusions 

The Mathematica package InvariantsSymmetries.m presented in this chapter can be 
used for computer-aided integrability detection of systems of nonlinear PDEs as they 
occur in various branches of science and engineering. 

More precisely, our package is a tool to search for the first half a dozen conservation 
laws and symmetries. If our programs succeed in finding a large set of independent 
conservation laws or symmetries, there is a good chance that the system has infinitely 
many of these quantities. For instance, if the number of conservation laws is 4 or 
less, most likely the system is not integrable — at least not in its current coordinate 
representation. 

Applied to a system with parameters, our package can determine the conditions 
on the parameters so that the system admits a sequence of conserved densities or 
generalized symmetries. 

An actual proof of integrability, by showing the existence of an infinity of 
conservations laws or symmetries, must be done analytically (see for results in 
this direction). On the other hand, constructing the recursion operator, and showing 
that it indeed satisfies the defining equation, provides conclusive proof of integrability. 
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Table 3 List of Software and Contact Information 


Name & System 


Scope 


Developer (s) & Address 


Email Address 


CONLAW 1/2/3 
(REDUCE) 


Conservation 
Laws 


T. Wolf et al. 
School of Math. Sci. 
Queen Mary & 
Westfield College 
University of London 
London El 4NS, U.K. 


T . Wolf (Smaths . qmw .ac.uk 


DELiA 
(Pascal) 


Conservation 
Laws and 
Generalized 
Symmetries 


A . Bocharov et al. 
Saltire Software 
P.O. Box 1565 
Beaverton, OR 97075 
U.S.A. 


alexeibQsaltire . com 


FS 

(REDUCE) 


Conservation 
Laws and 
Generalized 
Symmetries 


V. Gerdt & A. Zharkov 
Laboratory of Computing 
Techniques & Automation 
Joint Institute for 
Nuclear Research 
141980 Dubna, Russia 


gerdtO j inr . dubna . su 


Invariants 
Symmetries .m 
(Mathematica) 


Conservation 
Laws and 
Generalized 
Symmetries 


U. G6kta§ & W. Hereman 
Dept. of Math. Comp. Sci. 
Colorado School of Mines 
Golden, CO 80401, U.S.A. 


unalgSwolf ram . com 
wheremanOmines . edu 
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Table 3 cont. List of Software and Contact Information 


Name & System 


Scope 


Developer (s) & Address 


Email Address 


SYMCD 
(REDUCE) 


Conservation 
Laws and 
Generalized 
Symmetries 


M. Ito 

Dept . of Appl . Maths . 
Hiroshima University 
Higashi-Hiroshima 
724 Japan 


itoOpuramiis . arniath . 
hiroshima-u .ac.jp 


symmetry & 
mast er symmetry 
(MuPAD) 


Generalized 
Symmetries 


B. Fuchssteiner et al. 
Dept . of Mathematics 
Univ. of Paderborn 
D-33098 Paderborn 
Germany 


bennoSuni-paderborn . de 


Tests for 
Integrability 
(Maple & FORM) 


Conservation 
Laws , Genera- 
lized Symmetries, 
and Recursion 
Operators 


J . Sanders & J .P . Wang 

Dept. of Math. 

& Comp. Sci. 

Vrije Universiteit 

1081 HV Amsterdam 

The Netherlands 


jansaScs .vu.nl 


Tools for 
Conservation Laws 
(Maple) 


Conservation 
Laws 


M. Hickman 

Dept. of Maths. & Stats. 
University of Cainterbury 
Private Bag 4800 
Christchurch 
New Zealand 


M.Hickman 

(Smath . canterbury . ac . nz 
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